import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps

a = np.loadtxt('NWDOS_unconstraint_600K')
b = np.loadtxt('NWDOS_fix_PS4_600K')
c = np.loadtxt('NWDOS_fix_host_600K')
d = np.loadtxt('neck2_unconstrained.txt')
e = np.loadtxt('neck2_fixed_PS4.txt')
#print(simps(a[:,5],a[:,0]))

l = 6.59269


fig, ax1 = plt.subplots()
left,bottom,width,height = [0.58,0.22,0.3,0.3]
ax2 = fig.add_axes([left,bottom,width,height])

#ax1.title('Li NWDOS',fontsize=20)
ax1.plot(a[:,0]*2,a[:,1]/2,label='unconstrained',color='C1',zorder=5)
ax1.plot(b[:,0],b[:,1],label=r'fixed PS$_4^{3-}$',color='blue',zorder=1)
ax1.plot(c[:,0],c[:,1],label='fixed host',color='black',zorder=0)
ax1.set_xlim(0,40)
ax1.set_ylim(0,0.018)
ax1.set_yticks([0.00,0.009,0.018])
ax1.minorticks_on()
ax1.set_xlabel('Energy (meV)',fontsize=20)
ax1.set_ylabel('DOS (1/meV)',fontsize=20)
ax1.tick_params(which='both',direction='in',labelsize=20)
ax1.legend(frameon=False,fontsize=20,loc='upper left',handlelength=1,handletextpad=0.4)


ax2.plot(d[:,0],d[:,1],color='C1',zorder=5)
ax2.plot(e[:,0],e[:,1],color='blue',zorder=1)
ax2.axvline(l,0,0.5,color='black')

ax2.set_xlim(5,8)
ax2.set_ylim(0,2)
ax2.set_xticks([5,6,7,8])
ax2.minorticks_on()
ax2.set_xlabel(r'Area ($\rm{\AA}^2$)',fontsize=16,labelpad=-4)
ax2.set_ylabel('Bottleneck \n area dist. (a.u.)',fontsize=16,ha='center')
#ax2.set_ylabel('Bottleneck area distribution (a.u.)',fontsize=20)
ax2.tick_params(which='both',direction='in',labelsize=16,pad=4)

plt.savefig('Constraint_DOS_bottleneck.eps',format='eps',bbox_inches='tight')
plt.show()
